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ABSTRACT 


The IPS (Instantaneous Power Spectrum) spectral analysis technique has been the 
subject of study for many years. This thesis implemented the IPS algorithm using 
MATLAB. In addition, two additional programs were written to deal with progressively 
larger data sets. Based on a third order cumulant, the 1 % D spectral analysis technique, 


thought to perform well in low signal to noise environments, is also explored. 
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L STATIONARY SPECTRAL TECHNIQUES 


A. INTRODUCTION . 

In the disparate fields of geology, communications, astronomy, oceanography, 
chemistry and biomedicine, observations of physical processes are typically collected and 
then analyzed to extract as much information as possible. These observations are 
primarily analog and continuous in time. Spectral analysis is one of the initial analysis 
tools available; it is used to determine what frequency components are present in the 
observation. This analysis includes information about the spectral magnitude of the 
frequency components. Many spectral analysis techniques have been developed to 
investigate this problem. The first scientific application of spectral analysis, was in 
chemistry and astronomy. Light, whether from an astral object or from the combustion of 
a chemical sample was split into its component parts by prisms and the resulting spectra 
were photographed for later analysis [Ref. 1]. With the development of instruments which 
could record electromagnetic, seismic or acoustic information, additional spectral analysis 
techniques were developed to analyze these recorded observations. 

Classical spectral analysis techniques are based on Fourier methods and include the 
Periodogram [Ref. 2] and its time-frequency analog, the Spectrogram [Ref. 4], both of 
which will be discussed further in this chapter. The Fourier based methods are particularly 
well suited as computer-based analysis tools. The continuous-time data sequence is 
appropriately filtered and sampled to produce discrete data elements which can then be 
easily processed with the Fast Fourier Transform (FFT), a fast version of the Discrete 


Fourier Transform (DFT). 





Spectral analysis techniques based on linear filter models were also developed 
[Ref. 3]. The goal of these methods is to design a linear filter, H(e)), which when driven 
by stationary white noise can produce the data in a statistical sense. As an example, the 


Autoregressive (AR) model will return an all-pole filter with a resultant spectrum defined 





as 
2 
§ (ec) =—Pel__. (1.1) 
[A(e’*) 


where A(e’) is the Fourier transform of the pole coefficients and 5, represents the DC 

- gain of the filter. The AR model is well-suited to the detection of discrete sinusoidal 
signals or other narrow band signals. Other model-based spectral estimation techniques 
include both the Moving Average (MA) and the Autoregressive-Moving- Average 
(ARMA) models [Ref. 3]. Each of these has advantages and disadvantages depending on 
assumptions that can be made concerning the data sequence. All of these model-based 
methods, however, assume that the data sequence is Wide Sense Stationary (WSS). Many 
data sequences, however, are not WSS that is, the signals of interest are transient in a time 
domain sense or have dynamic spectral features. In these cases, the AR, ARMA and MA 
spectral techniques will not be adequate. 

Other spectral estimation techniques are based on eigenvalue decomposition or 
subspace methods (Ref. 7]. These methods presuppose that the data can be separated into 
a noise subspace and a signal subspace. The object, then, is to correctly identify which 
eigenvalues belong to the noise subspace and which belong to the signal subspace. 
Spectral techniques based on subspaces include Pisarenko, MUSIC and ESPRIT 
techniques [Ref. 3]. These techniques do not actually provide a classical spectral display, 


their spectra consists of delta functions representing the sinusoidal components. One 





major weakness of these methods is the fact that they were designed to isolate narrow 
band signals and are less suited for wide band signals. 

This thesis will explore two other spectral estimation techniques which are based on 
instantaneous estimates of the second and third moments. The first of these techniques to 
be explored is the Instantaneous Power Spectrum (IPS) [Ref. 4] which is based on an 
estimate of the second moment. It differs from the Periodogram, which is also based on 
an estimate of the second moment, by using an instantaneous estimate of the 
autocorrelation sequence. IPS belongs to a class of distributions called "Cohen's class" 
[Ref. 3]. The other technique to be explored is the 1 % D Instantaneous Power Spectrum 


(1 %D) [Ref. 6], which is based on third moment properties of the data. 


B. CLASSICAL SPECTRAL ESTIMATION TECHNIQUES 

All spectral estimation techniques are concerned with determining the spectral 
content of a finite set of observations. The formal definition of the Power spectral density 
(PSD) [Ref. 3] of a Wide Sense Stationary (WSS) discrete time random process is 


P(e!?)= Dre (Ae “TC WET; (1.2) 


k=~0 


where r,. (k) is the autocorrelation function of the observation, x(n), defined as 


re (k) = E[x"(n) x(n + 4) | (1.3) 


and where E denotes the expectation operator. To obtain the expected value, time 
averaging is usually used requiring essentially the observation of an ergodic process over 
an infinite period. The PSD displays the frequency components of any WSS random 


process of infinite duration. However, data observations are always of finite duration and 


wide-sense-stationary only in a local sense. The challenge then, is to approximate the true 








PSD with as much fidelity as possible. Classical spectral estimation techniques include the 
Periodogram and its time-frequency analog, the Spectrogram. 
1. Periodogram 


The Periodogram is defined as 


P(e) =|x(e”)) | (1.4) 


where P, (e’”) is the estimated PSD and X(e’®) is the Fourier transform of the 
observation x(n) of length Ny. One of the inherent weaknesses of the Periodogram is the 
fact that it deals with a finite set of observations. A finite set of observations can be 
obtained from an infinite set by applying a rectangular window of the form 


1, lsksN, . a5) 
0; otherwise 


nor| 


This time domain operation appears as a periodic convolution in the frequency domain. 
That is, the Fourier transform of the rectangular window is convolved with the Fourier 
transform of the observation set. The Fourier transform of the rectangular window is the 
digital sinc function which tends to smear the true PSD. To illustrate this point we can 
create an analytic sinusoid at some arbitrary frequency whose true Fourier transform 
would be an impulse function at the appropriate frequency location. For Figures 1.1 
through 1.3, an analytic signal in additive white Gaussian noise with a signal to noise ratio 
(SNR) of -10 dB was created. The signal to noise ratio is defined as 


Power, 
SNR = 1Orlog| Sera | oe) 


Power 


noise 


Lv) 


is 





Periodograms were then taken with the following parameters: . 
e Figure 1.1, 64 point observation sequence using a 64 point Fourier Transform 
e Figure 1.2, 64 point observation sequence using a 128 point Fourier Transform 
e Figure 1.3, 64 point observation sequence using a 512 point Fourier Transform 
Each Periodogram was plotted with a discrete plot function to enhance the display of the 


sinc function. In each figure the true frequency location of 15.78 is indicated by a straight 
line at the appropriate location 
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Figure 1.1. 64 point Periodogram 
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Figure 1.2. 128 point Periodogram 
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Figure 1,3. 512 point Periodogram 
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A comparison of Figures 1.1, 1.2 and 1.3 shows that as the transform length is 
increased, hence the number of frequency bins is increased, the position of the dominant 
spectral peak approaches more closely the true location of the frequency of the analytic 
signal. It can also be noted from these figures that although the apparent resolution 
improved as the transform length increased, the variance of the additive Gaussian white 
noise stayed at the same level. For the results in Figures 1.2 and 1.3, zeros were padded | 
to the data sequence to obtain the transform lengths of 128 and 512, respectively. The is on 
variance of the spectral estimate is independent of the length of the data sequence. 
Typically, sequential periodograms are taken and then averaged to reduce the undesirable 
effects of the variance of the estimate [Ref. 2]. Figure 1.4 is a 512 point Fourier transform : : 


Periodogram of a noise-free data sequence created with two analytic sinusoids separated 
by just oath of the sampling frequency. As in the other figures, the true frequency 


locations are indicated by straight lines. It can be seen that the Periodogram correctly 
resolves the two narrow-band components, but even in this ideal, noise-free environment, 


the position of each spectral peak is slightly off the true frequency bin locations. 
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Figure 1.4. 512 point Fourier Transform Periodogram of Two Analytic Sinusoids in additive 
Gaussian White Noise at a SNR of -10 dB 





Another weakness of the Periodogram is its inability to provide any 
information on the occurrence in time of an observation. To illustrate this an FSK data set 
composed of an analytic sinusoid which is switched from one frequency to another 


frequency at time 64 was created and is shown in Figure 1.5. 
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Ovpaervation Time 
Figure ].5. Time plot of Frequency Shift Key (FSK) signal 
When the Periodogram (shown in Figure 1.6) is calculated for this FSK observation, the 
presence of both sinusoids is evident but the fact that they existed during different times in 


the observation set cannot be seen. 
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Figure 1.6, 128 point Fourier Trancfoen Pericdograil of a FSK Signal 
The ameaequeney representation of the Periodogram, the Spectrogram, can show that 
these signals were separated in time. 
2. Spectrogram 
The Spectrogram operates by applying a window to a subset of the data 
set and computing and saving the Periodogram of that subset. Then the window is 
stepped through the data by some fixed numbers of points and a Periodogram is 
recomputed. This process of stepping through the data and sequential computation of 
periodograms is repeated until the desired number of spectrogram lines is obtained. 
Figure 1.7 is a mesh and contour plot of the spectrogram of the FSK data set of Figure 
1.5. 

















(a) p 4 @) 
Figure 1.7. Mesh and Contour plots of the Spectrogram of a FSK Signal 
It can be seen that the frequency location is most concentrated at the 
beginning and end of the surface. At these points, the window includes a data segment 
which contains only one of the signals. As the window is stepped through the data set, the | 
number of data points that include one frequency or the other changed. The Spectrogram 
surface seems to indicate that the first frequency at frequency location 5 is present from 
time 0 to time 41 when in fact it was only present from time 0 to 32. On the other hand, 
the second signal at frequency location 20 begins at time 32 not at time 26 which seems to 
be indicated. | 
~ An additional test signal whose frequency changes linearly with time, 
such as the linear FM Chirp, was created. Figure 1.8 is the Periodogram of a linear FM 
Chirp. 
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Figure 1.8. Periodogram of a Linear FM Chirp 
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The Periodogram of Figure 1.8 indicates the presence of several sinusoidal components 
rather than the presence of a linear FM chirp. The Spectrogram in Figure 1.9, however 





r 
clearly shows that the data is a linear FM chirp 
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(a) a 
Figure 1.9. Spectrogram of a Linear FM Chirp 





Using the linear FM chirp as an example ofa dynamic signal, it can be seen that the width — 
of the spectral lobe is wider than the spectral lobes for the FSK signal in Figure 1.7. To 
further illustrate this phenomenon a data set composed of a signal whose frequency 
changes as a quadratic function of time, a quadratic FM chirp, from frequency location 1! 


to location 30 was created. Figure 1.10 is the spectrogram of this signal. 

























A 
ry 






ee ar 





(a) | 7s | (b) 
Figure 1.10. Spectrogram of a Quadratic FM Chirp 
The true position of the instantaneous frequency can be seen as a hyperbolic curve on the 
contour subplot of Figure 1.10(b). The quadratic FM Chirp is an even more dynamic 
signal than the linear FM chirp. The spectral lobe broadens, best seen in Figure 1.10 (b), 
as the signal's frequency change accelerates. 
We want, therefore, to look at other spectral techniques which may 


moderate some of the inherent weaknesses of the Spectrogram. 


i. INSTANTANEOUS POWER SPECTRUM 


A. INTRODUCTION 
The Instantaneous Power Spectrum (IPS) is based on a class of time-frequency 
distributions called "Cohen's class" [Ref. 3]. Cohen's class is defined for a continuous 


signal by 


c=] ffowne+ Ser e-Demeavatae 20 


Many different power spectral techniques can be derived from this class of distributions 
including Wigner-Ville, IPS and Rihaczek time-frequency representations [Ref. 4]. Each 
of these techniques can be obtained from the generalized expression by the selection of an 
appropriate kernel, ¢.(v,t). In the case of IPS, the kernel used is e’””. Originally the 
estimate of the instantaneous frequency content utilizing Page's definition of the 


instantaneous power spectrum as the derivative of a running spectrum [Ref. 8]: 


x | Fe EPR: 
PUS=S|S, 1) BD 
where 


S>(f)= [s(cjeP* ar. (2.3) 


The concept of the instantaneous power spectrum was further expanded by Levin [Ref. 9] 


with the addition of a backward running spectrum defined as: 





PUN= Zh") (2.4) 


where 


S(f)= | s(t)eP "dr. | (2.5) 
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Using both the forward and backward running spectrums we can define IPS as the average 


of these two spectra given by 


PSA\=S[P UA UA) (2.6) 


The discrete form of the IPS algorithm utilized in this thesis is [Ref. 5] 


‘  . 
IPS_(n,00) == Ye {x(n)x" (n= k) +x" (n)x(n+ wow (ye. (2.7) 


where w(k) is a window function, x(7) is the data sequence and N., is the length of the 
data sequence. The discrete IPS expression is actually the Discrete Founer Transform 


(DFT) of 


{x(n)x"(n— k)+x"(n)(n+k)}po(O)w(k)_ (2.8) 


In order to exploit the efficiency of the FFT the length of the data sequence is constrained 
to be a power of 2, such as 64, 256, or 512. 

This chapter will explore three types of IPS programs which are given in the 
Appendix. The IPS program is designed for the analysis of relatively short data 
sequences, i.e., 64 to 1024 points. The IPSSURF program should be selected when ee 
lengths are between 1024 and 213 to 215 points. The IPSLOFAR program should be 


selected to process data sets larger than 215 points. 


B. TRADITIONAL IPS TIME-FREQUENCY DISPLAY 

For short data records, the traditional time-frequency display of the IPS algorithm is 
well-suited. Input parameters of the program define the dimensions of the returned 
time-frequency surface. The parameters include the window type (Rectangular or 


Hamming), the window length (normally half of the data sequence length) and the step 


14 


ir 








(the distance through which the window is stepped through the data sequence). As an 
example, if the data sequence is 512 points long with a window length of 256 points and a 
step of 1, the resulting surface contains 512 rows in the time direction having a frequency 
range of -% to +n divided into 256 frequency bins. The program assumes the data 
sequence to be of fixed duration. Since the IPS algorithm looks both backward and 
forward in time, the data sequence is padded with zeros at the beginning and end equal to 
the width of the selected window. The program returns only the positive frequency half of 
the resulting time-frequency surface to limit the size of the final display. In Figures 2.1 
through 2.6, test data sets are used as signals of length 128. The window length is chosen 
to be 64 points with a Hamming window and the step is chosen to be 1. 

1, Single Analytic Sinusoid 


A single analytic sinusoid was created whose frequency location is exactly 19, 
‘ nee 19 ee: . 
i.e., the digital frequency was rae This implies that for a window length of 64, the 


spectral peak should occur at bin 19. Had the window length been 128 points, the 


frequency location would have been 38. The IPS surface, both mesh and contour 


subplots, are shown in Figure 2.1. 
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Figure 2.1. Single Analytic Sinusoid via IPS 
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It can be scenshek the frequency location is symmetrically centered at the proper location 
of 19. The width of the spectral ridge in Figure 2.1(a), measured at the point where the 
magnitude has dropped by 3 dB, is approximately three frequency bins wide. For this 
single analytic sinusoid data set, the sensitivity of the IPS algorithm is comparable to the 
Spectrogram. The IPS builds more quickly to a constant amplitude and trails off less 
quickly at the end of the data set than the Spectrogram would for the same data set. The 
sidelobes, evident on the mesh subplot, are never of such amplitude that they obscure the 
position of the spectral ridge and have dampened completely as the window moves well 


into the data set. 


2. Multiple Analytic Sinusoids 
A data set consisting of two analytic sinusoids was created with frequency 
locations of 19 and 25, respectively. The IPS surface, both mesh and contour subplots is 


shown in Figure 2.2. 








(a) a  @) 


Figure 2.2. Multiple Analytic Sinusoid via IPS 
It can be seen that the spectral ridges are symmetrically centered at the proper locations of 


19 and 25 and their 3 dB widths are, as in Figure 2.1, still approximately three frequency 


16 





bins wide. The sidelobes, evident on the mesh surface, dampen as the window igves into 
the data set. The modulation, evident along the spectral ridges, is a consequence of 
cross-spectral terms which ride on the autocorrelation terms. One significant advantage of 
the IPS algorithm, compared to the Wigner-Ville algorithm, is that it does not experience 
spurious cross-spectral terms between the strong spectral peaks. The intra-ridge 
modulation of the IPS affects the display of the time-frequency surface but does not 
degrade or obscure the determination of the spectral locations. 

3. Frequency Shift Key (FSK) 

A Frequency Shift Key (FSK) data set was created composed of an analytic 


. : : 10 : 
sinusoid whose frequency was shifted from a to a midway through the data set of 128 


points. The resultant IPS surface, both mesh and contour subplots is shown in Figure 2.3, 


the corresponding Spectrogram is shown in Figure 2.4. 








(a) | 6) 


Figure 2.3. IPS of FSK data set 
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(a) (b) 
Figure 2.4. Sicciontain of FSK dataset | | 
Both the IPS surface and the Spectrogram correctly locate the respective frequency 
locations, however the IPS surface accurately discerns the time of the signals. The IPS 
surface clearly shows that the frequency shifted at time 64 to a higher frequency. The 
Spectrogram incorrectly shows an overlap in the time of the two signals. The IPS surface 
does éxperience sidelobes which never completely disappear through its extent, however it 


can be argued that the correct spectral location in frequency is of more significance. 
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4. Linear FM Chirp 


A data set consisting of a signal whose frequency changes linearly with time, 


i.e. linear FM chirp, was made to transition from a frequency of a to ~ The resultant 


IPS surface, both mesh and contour subplots, is shown in Figure 2.5. 





. 7 Figure 2.5. IPS of Linear FM Chirp 
The true position of the instantaneous frequency can be seen as a straight line on the 
contour subplot of Figure 2.5(b). The width of the spectral ridge is broader than the 
previous figures, however, drawing a line along the centerline of the structure would 
clearly correctly locate the instantaneous frequency at the correct time in the data set. 
5. Quadratic FM Chirp 
A data set consisting of a signal whose frequency changes as a quadratic 


function of time, i.e. over the length of the data the quadratic FM chirp, was made to 


0 
transition from a frequency of =, to = The resultant IPS surface, both mesh and 


contour subplots is shown in Figure 2.6. 














each of the sinalter segments and then concatenating the surfaces together into a larger 
surface. Recall that the IPS program pads the data sequence with zeros to allow the 
algorithm to move backward and forward in time. The IPSSURF program does pad the 
first and last smaller data segments with zeros but the other data segments are padded 
with true past and future data points from the full data sequence. The IPSSURF 
parameters include the window type (Rectangular or Hamming), siglen (the desired length 
of the smaller data segments), the window length (normally half of the smaller data 
segment length) and the step (the distance through which the window is moved through 
the smaller data segments). As an example, if the original data sequence is 2048 points 
long and smaller data segments of 512 were selected with a window length of 256 points 
and a step of 8, the resulting surface contains 256 rows in the time direction having a 
frequency range of -1 to + divided into 256 frequency bins. The IPSSURF program, if 
invoked with a step size of 1, would concatenate full IPS surfaces into a very large 
surface. For this reason, the IPSSURF program is designed to be used as a broader 
analysis tool for an overall look of a large data sequence by invoking it with a step size 
larger than 1. Once an area of interest has been isolated with the IPSSURF program, finer 
analysis would be done using the IPS program. The test signals created for Figures 2.8 to 
2.11, were 1,024 points long, the smaller data segment length (the siglen parameter) was 
selected as 128 points, the window length was selected as 64 points and the step was 
selected as 8. The IPSSURF subplots for FSK signal (Figure 2.8); Linear FM chirp 
(Figure 2.9); Quadratic FM chirp (Figure 2.10) and the Multiple Component Signal 
(Figure 2.11) follow. Each of the data sets were created as defined in sections B.3 


through B.6. - wih gs, ha oe 
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Figure 2.11. Multi-component Signal with IPSSURF 

Comparison of Figures 2.8 through 2.11 with the corresponding IPS Figures 2.3, 
2.5 through 2.7 seem to indicate that the IPSSURF program shows finer spectral ridge 
details than the IPS surfaces. It can also be seen that the sidelobes on the IPSSURF 
surfaces seem to dampen more quickly than the IPS surfaces. These effects, however, are 
a result of time compression, determined by the selection of the step size parameter. As 
the step size is increased, finer details of the spectral surface would be progressively 
obscured. For this reason, the IPSSURF program should be used as a coarse analysis tool 
to locate and define areas of a large data set which could be further analyzed with the IPS 


program. 


D. IPSLOFAR . 

The IPSLOFAR program (Appendix) is for very large data sets. IPSLOFAR is 
based on a "waterfall" display routinely used in the display of sonar data, called the 
LOFARGRAM. The IPSLOFAR program, like the IPSSURF program, divides a large — 
data set into smaller segments and calculates IPS surfaces for each of the smaller 


segments. Unlike the IPSSURF program, the average over time is then calculated and 
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placed as a row ina larger surface for display. The IPSLOFAR program can be seen as a 
part of the continuum of ‘IPS’ programs where the IPS program is the finest analysis tool 
and the IPSLOFAR program is the coarsest analysis tool. The IPSLOFAR program 
parameters include the window type (Rectangular or Hamming), siglen (the desired length 
of the smaller data segments), the window length (normally half of the smaller data 
segment length) and the step (the distance through which the window is stepped through 
the smaller data segment). 

A data set was created of Jength 16,384 points. The signal is composed of a linearly 
decreasing FM chirp made to transition from frequency location 30 to 1, a two component 
FSK signal which switched from frequency location 10 to 20 midway through the data 
sequence and an analytic sinusoid at frequency location 30. The IPSLOFAR parameters 
used were a siglen of 128 points, a window length of 64 points, a Hamming window and a 
step of 8. The IPSLOFAR contour plot is shown in Figure 2.12(a) and a plot of the 
average over time of the surface is shown in Figure 2.12(b). 

The IPSLOFAR display clearly shows all signal components. The width of the 
spectral ridges is comparable to both the IPS and IPSSURF programs. The Lofargram, 
used in sonar displays, is actually an intensity plot, vice a contour plot. Because we can 
only use contour plots for IPSLOFAR in Matlab, the display is somewhat different from a 


traditional Lofargram. 
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Figure 2.12. IPSLOFAR surface for a Multi-component Signal 
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HL 1% D INSTANTANEOUS POWER SPECTRUM 


A. INTRODUCTION 

Both the Spectrogram and IPS techniques are second moment spectral analysis 
techniques. They and other related techniques are well-suited to extract spectral 
information, however when additive Gaussian noise perturbs signals, the application of 
higher order moments is thought to improve the signal to noise performance of the 
spectral algorithm. For zero mean Gaussian random processes, the third moment and 
third order cumulant are equal to zero. The 1 % D instantaneous power spectrum 
[Ref. 6] is based on a third order cumulant. Cumulants are related to ordinary moments, 
as a matter of fact, for a zero mean random process the first and second order cumulants 
are equal to the first and second moments. For higher order spectral techniques, the use 
of cumulants is often preferred to higher order moments because cumulants can measure 
the departure of a random process from a Gaussian random process [Ref. 3]. The 1 4D 
instantaneous power spectrum implemented in this thesis is derived from the bispectrum. 
The bispectrum is a third order spectrum defined as [Ref. 3] 


SOM e)= Ye VCO TU be Orn; (3.1) 


h=~t I= 
where C,°’[],, 1, ] is the third order cumulant. If we assume that we have a zero-mean 


random process the third order cumulant is 


COTA] =E{x'(n)x(n+),)x(nt1)}. (3.2) 


We can derive a degenerate estimate of this cumulant by setting /; to zero, therefore 


C,° [0,1] = E{x*(n)x(n)x(n+1,)} | (3.3) 
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and 
C0, = EXlx(f x(n +}. (3.4) 
Furthermore, when we replace the expectation with the instantaneous value the third order 
spectrum becomes 


S,(e’”) = ¥ [xn [x(n ye (3.5) 


kane : 
As in the derivation of the IPS technique of Chapter II, we will utilize both Page's | 
definition of the instantaneous power spectrum and Levin's addition of a backward running 
spectrum to finally derive the discrete form of the 1 % D instantaneous power spectrum 
used in this thesis as 
N, 


1% Do) = + lletol (0a) +e xtnt k) hw (O)w(ke/ - (3.6) 


where w(k) is a window function, x(7) is the data sequence and N, is the length of the 
data sequence. As in the IPS technique, the length of the data sequence is constrained to 
be a power of 2, such as 64, 256, or 512. | 

The three 1 % D programs are given in the Appendix. The ONE_HALF program is 
used for relatively short data sequences, i.e., 64 to 1024 points. The ONESURF program 
is used for data sequences typically between 1024 and 2!3 and 215 points. The 
ONELOFAR program should be used for data sets larger than 2 13 points. The 
philosophy of the 1 % D programs is analogous to the corresponding IPS programs of 


Chapter IT. 


B. TRADITIONAL 1 % D TIME-FREQUENCY DISPLAY 
The parameters for the ONE_HALF program include the window type (Rectangular 


or Hamming), the window length (normally half of the data sequence length) and the step 
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(the distance through which the window is stepped through the data sequence). The 

program returns only the positive frequency half of the resulting time-frequency surface to 
limit the size of the final display. For Figures 3.1 through 3.6, test data sets of length 128 
were created. The window length is chosen to be 64 points with a Hamming window and 
the step is chosen to be 1. Figures 3.1 through 3.6 utilize the same data sets as were used 


in Figures 2.1 through 2.3 and 2.5 through 2.7 of Chapter II, respectively. 


1. Single Analytic Sinusoid 
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(a) | iG (b) 
Figure 3.1. Single Analytic Sinusoid via 1% D 

When compared to Figure 2.1, the appearance of the spectral ridge on the 
mesh subplot of Figure 3.1, is not as smooth in appearance. The width of the spectral 
ridge measured at the point where the magnitude has dropped by 3 dB is comparable to 
the IPS technique. The sidelobes on the time-frequency surface are less defined but the 
extent of their effect is again comparable to the IPS technique. It can also be seen that the 
1 % D technique quickly reaches a constant amplitude and falls off in times comparable to 
the corresponding IPS case. Also during the central region of the surface, along the time 
axis, the noise background seems to be smaller in Figure 3.1 than the corresponding 


Figure 2.1. 


29 








2. Multiple Analytic Sinusoids — 











(a) ee (b) 
Figure 3.2. Multiple Analytic Sinusoid via 1 4D | 
When compared to the corresponding IPS surface, Figure 2.2, the 1 4D 
time-frequency peaks at the true frequency locations having comparable spectral ridge 
widths. The envelope of the spectral ridges fluctuates more than on the corresponding 
IPS surface. 
3. Frequency Shift Key (FSK) 





eso H 






in Cee he 





Time 


£ h 

its 
PPS 
Fina 





Frequency Bie 


(a) . a (b) 
Figure 3.3. FSK signal via ONE_HALF 
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For the FSK signal data set, areas in the 1 % D surface where the frequency 


shift occurs differs from the corresponding IPS surface of Figure 2.3. The 14D 


technique does not delineate the time of the signals as accurately as the IPS technique. 


The time of the frequency shift is obscured by the appearance of cross terms. 


4. Linear FM Chirp 
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(a) 


Figure 3.4. IPS of Linear FM Chirp 
The true position of the instantaneous frequency can be seen as a straight line 
on the contour subplot of Figure 3.4(b). The width of the 1 % D spectral ridge is wider 
than the corresponding IPS surface in Figure 2.5. The surface has a coarser appearance, 
showing less detail, than the corresponding IPS surface. The envelope of the spectral 
ridge shows more variability with slower intra-ridge modulation (1.e., fluctuations along 


the ridge top are at a lower frequency) than the corresponding IPS surface. 
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5. Quadratic FM Chirp 
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Figure 3.5. IPS of a Quadratic FM Chirp 

The true position of the instantaneous frequency is indicated by a hyperbolic 
line on the subplot of Figure 3.5(b). Once again, when compared with the IPS surface in 
Figure 2.6, it can be seen that the width of the spectral ridge is broader than in the IPS 
case, however the true instantaneous frequency can be obtained as an average along the 
spectral ridge. It can also be seen that the amplitude of the spectral ridge decreases 
towards the end of the data set and also broadens. This example shows that the 1 % D is 
less suited for a dynamic signal than IPS but better suited than the spectrogram, if one 


disregards issues such as noise. 
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6. Analytic Linearly Increasing FM Chirp, Quadratically Decreasing FM 


= Chirp and stationary Sinusoid (Multiple Component Data Set) 
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(a) = 4 (b) 
Figure 3.6. IPS of Multiple Component Signal 
When compared with the IPS surface in Figure 2.7, the width of the spectral 
ridges indicating the dynamic signal components, i.e., the linear FM chirp and the 
quadratic FM chirp is wider than IPS. The width of the spectral ridge which indicates the 
stationary sinusoid is comparable to the IPS surface, although there is more modulation 


apparent on the spectral ridge. 


C. LINKED 1 % D SURFACES 

The parameters for the ONESURF sroavim include the window type (Rectangular 
or Hamming), siglen (the desired length of the smaller data segments), the window length 
(normally half of the smaller data segment length) and the step (the distance through 
which the window is moved through the smaller data segments). The ONESURF program 
is designed to be used as a broader analysis tool for an overall look of a large data 
sequence by invoking it with a step size larger than 1. The ONESURF program is used to 
plot the 1 % D surface of an FSK signal (Figure 3.7); a linear FM chirp (Figure 3.8), a 
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quadratic FM chirp (Figure 3.9) anda multiple component signal. Each of the data sets 


were created as defined in Chapter II, sections B.3 through B.6. 



































(a) (b) 
Figure 3.7. FSK via ONESURF 
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Figure 3.8. Linear Chirp via ONESURF 
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(a) 
Figure 3.10. Multi-component Signal with IPSSURF 


Comparison of Figures 3.7 through 3.10 with the corresponding IPSSURF surfaces 
in Chapter II, Figures 2.8 through 2.11 shows very few differences. The ONESURF 
program, as the IPSSURF program, is best suited as a coarse analysis tool to locate and 
define areas of a large data set which could be further analyzed with the ONE-HALF 


program. 
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D. ONELOFAR 
The ONELOFAR program parameters include the window type (Rectangular or | ; 
Hamming), siglen (the desired length of the smaller data segments), the window length | af 
(normally half of the smaller data segment length) and the step (the distance through 
which the window is stepped through the smaller data segment). 
The experimental test signal used is the same as used in the IPS chapter, section D. 
The ONELOFAR parameters used are a siglen of 128 points, a window length of 64 
points, a hamming window and a step of 8. The ONELOFAR contour plot and 
corresponding plot of the average over time of the surface is shown in Figure 3.11(a) and 
(b), respectively. The ONELOFAR display clearly shows all signal components. The plot 
of the average over time also helps to locate the spectral components of the signal, 


especially in low SNR environments, That is, it can be used to extract information at very 


low SNR's, provided the spectral components are stationary. 
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Figure 3.11, ONELOFAR surface for a Multi-component Signal 
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IV. PERFORMANCE AND COMPARISON OF SPECTROGRAM, 
INSTANTANEOUS POWER SPECTRUM AND I % D TECHNIQUES 


A. INTRODUCTION 

Chapters I, II and III discussed, respectively, the Spectrogram, IPS and 1'%D 
techniques. Several test data sets were created and subsequently analyzed. Now we 
investigate the strengths and the weaknesses of each of these techniques in relation to one 
another with signals embedded in Gaussian white noise. A typical ocean acoustic data set 


is also processed. — 


B. SPECTRAL SENSITIVITY WITH SINUSOIDS IN ADDITIVE GAUSSIAN 
WHITE NOISE 
A data set was created, composed of sinusoids at frequency locations 5, 10, 15, 20, 
25 and 30 which were embedded, respectively, in additive Gaussian white noise at -6, -3, 
0, 3, 6 and 9 dB SNR. Figures 4.1 through 4.3 show the respective Spectrogram, IPS and 
ONE_HALF (1 % D) time-frequency surfaces in subplots using 
e (a) mesh surface 
e (b) contour plot with the true frequency locations shown by a line at each location 
e (c) plot of the average over time of the corresponding time-frequency surface. 
The integration time for the spectrogram was 64 while the window length for the 
IPS and 1 % D programs was also 64. The mesh subplots for both the IPS and 14D 
techniques show their characteristic modulation along the spectral ridges due to the 
cross-spectral terms. The Spectrogram mesh surface has no modulation along its spectral 
ridge, which is also characteristic for this technique. The respective contour plots of the 


three techniques offer more information. It can be seen that the width of the spectral ridge 
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is narrower in the Spectrogram compared to the other two techniques. However, it can 
also be seen that where the Spectrogram can locate the sinusoid in the Gaussian white 
noise to down to 3 dB SNR, the IPS techniques can locate down to -3 dB SNR and the 
1 4 D technique even locates to -6 dB SNR. 


39 















VSS 
‘ZS: 
SS 
=< 
————— 

















oO 


Figure 4.1. Spectrogram of Sinusoids embedded in Gaussian White Noise 
at -6, -3, 0, 3, 6 and 9 dB SNR respectively at frequency locations 
5, 10, 15, 20, 25 and 30 
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(c) 
Figure 4.2. IPS of Sinusoids embedded in Gaussian White Noise 


at -6, -3, 0, 3, 6 and 9 dB SNR respectively at frequency locations 
5, 10, 15, 20, 25 and 30 
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(c) 
Figure 4.3. 1 4 D of Sinusoids embedded in Gaussian White Noise 
at -6, -3, 0, 3, 6 and 9 dB SNR respectively at frequency locations 
5, 10, 15, 20, 25 and 30 
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C. SPECTRAL SENSITIVITY WITH AN OCEAN ACOUSTIC DATA SET 
Figures 4.4 through 4.6 show the respective Spectrogram, IPS and ONE_HALF 
(1 %D) time-frequency surfaces in subplots using 

@ (a) contour plot with the true frequency locations shown by a line at each location 

e (b) plot of the average over time of the corresponding time-frequency surface. 

The data set seems to be composed of several stationary signals embedded in a noisy 
background. As in the previous section, the Spectrogram does have narrower spectral 
ridges than both the IPSLOFAR and ONELOFAR time-frequency surfaces. On closer 
examination of the surfaces, note the spectral ridge at approximately frequency location 35 
in Figure 4.6. The ridge is also seen in the IPSLOFAR surface but is missing altogether in 
Figure 4.4, the Spectrogram surface. The 1 2 D surface does clearly locate a frequency 
component at location 35, but apparently at the expense of obscuring stronger signals in a 


very ‘busy’ surface. 
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Figure 4.5. Acoustic Data via PSLOFAR 
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Figure 4.6. Acoustic Data via ONELOFAR 
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D. CONCLUSIONS AND SUMMARY 

The Spectrogram is well-suited to quickly and clearly resolve stationary signals in 
moderately noisy environments. The IPS technique is superior to the Spectrogram both in 
the detection of stationary signals in noise and particularly in the detection of dynamic 
signals i.e., the Quadratic FM Chirp. Both the IPS and 1 % D techniques cannot resolve 
stationary spectral components as finely as the Spectrogram, but it could be argued that 
their superior ability to locate spectral components in noisier environments is more 
valuable in many cases. Potentially, one could also increase the integration time of the IPS 
based techniques to obtain better resolution. The judicious combination of these 
techniques can optimize the process of spectral analysis. The Spectrogram can be used as 
an efficient, first look at a particular data set. Areas of interest, which may indicate 
dynamic signals could be further analyzed by the IPS technique and finally the 1% D 


technique can offer superior location of spectral components in noisier environments. 


E. AREAS FOR FURTHER RESEARCH 

The derivation for the 1 4 D spectral technique used in this thesis may be further 
explored to extend its sensitivity to spectral components embedded in a noisy background. 
The bispectrum surface could be exploited by other than the 1 % D technique, which uses 
a degenerate form of the tricorrelation function. Applications of the Radon transform 
might shed new insights into the detection and classification of stationary and transient 
spectral components. Interpretation of each of these techniques depends on the display 
techniques which are available. The lofar-type displays for the IPSLOFAR and 


ONELOFAR programs could be enhanced if a true intensity plot could be obtained. 
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APPENDIX 
PROGRAM LISTINGS 
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IPS.M 


%function [P,freqindex,timeindex]=ips(data, wintype, winlen, step); 

% This function will calculate an Instantaneous Power Spectral (IPS) surface. 
%The IPS surface (output matrix P) characteristics are determined by the selection 
%of window type (wintype), window length (winlen) and the distance that the 
%window is moved through the data sequence (step). 

%The program plots only the positive half of the spectral plane. The 

% outputs timeindex and freqindex are the appropriate plot indices for the 
%output time-frequency surface. 


% 

%The inputs are: 

%data: The input data string 

%wintype: 'O' Rectangular Window 

% '!' Hamming Window 

%winlen: | The desired width of the window, normally half of the siglen 
%step: Time step desired, normally '1' or a multiple of '2' 

% 

%The outputs are: 

%P The IPS time-frequency surface 


%timeindex The y axis index 
Yfreqindex The x axisindex 
%See also IPSSURF, IPSLOFAR 


%Karen A. Hagerman 
%06 May 1992 


function [P,freqindex, timeindex]=ips(data, wintype, winlen, step) 
{m,n]=size(data); 


if m~=1 
data=data.'; 
end 


siglen=length(data); 
if wintype==0 
win=ones(winlen-1,1); 
elseif wintype==1 
win=hamming(winlen-1); 
end 
W=[win(winlen/2:-1:1)); 
x=[zeros(1,winlen) data zeros(1:winlen)]; 
p=zeros(siglen/step, winlen); 
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index=1; 

for n=winlen+1:step:siglen+winlen-stept+1 
Xm=[conj((x(n:-1:n-(winlen/2-1)))).’ (x(n:n+(winlen/2- 1))). 1; 
Xn=[x(n);conj(x(n))]; 
product=((Xm*Xn).*W).’; 
product=[product 0 conj(product(winlen/2:-1:2))]; 
p(index, :)=fftshift(real(.5*#R(product))); 
index=index+1, 

end 


p=p(:,winlen/2+1-:winlen); 
[prow,pcolumn]=size(p); 


%Smoothing 
p_temp(1,:)=mean(p(1:3,:)); 
p_temp(2,:)=mean(p(1:4,:)); 
p_temp(prow-1,:)=mean(p(prow-3:prow,:)); 
p_temp(prow, :}=mean(p(prow-2:prow,:)); 
for m=3:prow-2 

p_temp(m, i a m+2,:)); 
end 


pease 


freqindex=[0:pcolumn- 1]; 
timeindex=[1:prow]; 
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IPSSURF.M 


%function [Q,xindex, yindex]=ipssurf(sig, siglen, wintype, winlen, step); 

%This function will calculate an Instantaneous Power Spectral (IPS) surface made of 
smaller IPS surfaces. The smaller IPS surfaces are calculated for data sequence pieces 
%of length (siglen) of the total input data (sig). Characteristics of the smaller surfaces 
%are determined by the selection of window type (wintype), window length (winlen) and 
%the distance that the window is moved through the data sequence (step). The surfaces 
%are then appended edge to edge to form the output Q surface. 


% 
%The inputs are: 
Ysig: The input data string 


“%siglen: The desired length of the discrete pieces of sig 
Y%wintype: 'O' Rectangular Window 


% ‘l' Hamming Window 

%winlen: The desired width of the window, normally half of the siglen 
%step: Time step desired, normally '1' or a multiple of '2' 

% 

%The outputs are: 

%P The IPSSURF time-frequency surface 


%yindex = The y axis index 
Y%xindex The x axis index 
% 
%See also IPS, PSLOFAR 


%Karen A. Hagerman 
%06 May 1992 


function [Q,xindex, yindex]=ipssurf(data, siglen, wintype, winlen, step) 


if nargin==1 ; 
siglen=input(Enter the desired length of the sequence ‘); ; 
wintype=input(ENTER “0" for RECT. WINDOW or "1" for HAMMING WINDOW); 
winlen=input(‘ENTER length of the window (must be an even number):’); 
step=input(‘Input desired step in time '); 

end 


{m,n]=size(data); 
if m~=1 
data=data.'; _ 


end 


rows=floor(length(data)/siglen); 
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data=[zeros(1,winlen) data zeros(1,winlen)]; 
len=length(data), 

finish1=winlen/2-1, 

finish2=finish1+1,; 


if wintype==0 
win=ones(winlen- 1,1); 
elseif wintype==1 
win=hamming(winlen-1); 
end 


W=([win(winlen/2:-1:1)];; 
Q=zeros(rows*(siglen/step), winlen/2); 
qindex=1 :siglen/step:(rows* siglen/step); 
l=1; 
for k=1:siglen:len-siglen-2*winlen+1 
x=data(1,k:k+siglen+2*winlen-1); 
index=1; 
for n=winlen+1 :step:siglen+winlen-step+1 
Xm=[conj((x(n:-1:n-finish]))); (x(n:n+finish1))); 
Xn=[x(n) conj(x(n))]; 
product=((Xn* Xm).*W); 
product=[product 0 conj(product(finish2:-1: »)) 
P(index, :)=(fftshift(real(.5*fR(product)))); 
index=index+1; 
end 
Q(qindex(1):qindex(l)+(siglen/step)-1,:)=P(:, winlen/2+1:winlen), 
1=1+1; 
end 


Q=flipud(Q); 
[qrow,qcolumn]}=size(Q); 
xindex=[0:qcolumn-1]; 
yindex=[1:qrow]; 
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IPSLOFAR.M 


%function[Q, freqindex, timeindex]=ipslofar(data, siglen, wintype, winlen, step); 
%This function will calculate a 'Lofar' display for a selected data sequence. 
%The total data sequence is first divided into equal length pieces of 
%length (siglen). An IPS surface is then calculated for each piece. The 
%IPS surface characteristics are determined by the selection of window 
Ytype (wintype), window length (winlen) and the distance that the window 
%is moved through the data sequence pieces (step). The mean is then taken 
%of the IPS surface and is placed sequentially in the Q matrix for display. 
%The Q matrix plots only the positive half of the spectral plane. The 
“outputs timeindex and freqindex can be used in plots to interpret the 
“%results. 


% 
%The inputs are: 
Ydata: The input data string 


%siglen: The desired length of the discrete parts of the sequence 
Ywintype: 'O' Rectangular Window 


% ‘l' Hamming Window . 
%winlen: The desired width of the window, normally half of the siglen 
%step: Time step desired, normally '1' or a multiple of '2' 

% | 

“The outputs are: 

%P The IPSLOFAR time-frequency surface 


“timeindex The y axis index 
“ofreqindex The x axis index 
% 

%See also IPS, IPSSURF 


%Karen A. Hagerman 
%06 May 1992 


function (Q,freqindex,timeindex]=ipslofar(data, siglen, wintype, winlen, step) 


if nargin==1 
siglen=input(Enter the desired length of the sequence ‘); 
wintype=input(ENTER "0" for RECT. WINDOW or "1" for HAMMING WINDOW:’); 
winlen=input(ENTER length of the window (must be an even number): '); 
step=input(‘Input desired step in time '); 

end 

[m,n]=size(data); 

if m~=1 
data=data.'; 

end 

rows=floor(length(data)/siglen); 
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Q=zeros(rows, winlen/2); 

data=[zeros(1,winlen) data zeros(1,winlen)]; 

len=length(data); 

finish 1=winlen/2-1; 

finish2=finish1+1, 

if wintype== 
win=ones(winlen-1,1); 
elseif wintype==1 
win=hamming(winlen-1); 

end 

W=[win(winlen/2:-1:1)]'; 


qindex=1; 
for k=1:siglen:len-siglen-2*winlen+1 
=data(1,k:k+siglent+2* winlen-1); 
index=1; 
for n=winlen+1 :step:siglen+winlen-step+1 
Xm=([conj((x(n:-1:n-finish1))); (x(n: eens 
Xn=[x(n) conj(x(n))]; 
product=((Xn* Xm).*W); 
product=[product 0 conj(product(finish2:-1:2))]; 
P(index,:)=(fftshift(real(.5*fft(product)))); 
index=index+1; 
end 
if index ~=2 
Q(qindex, :=mean(P(:,winlen/2+1 :winlen)), 
qindex=qindex+1; 
else 
Q(qindex, :)=P(1,winlen/2+1: ume 
gindex=qindex+1; 
end 
end 


Q=flipud(Q); 
[qrow,qcolumn)=size(Q); 


freqindex=[0:qcolumn-1]; 
timeindex=[1:qrow]; 
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ONE_HALF 


%function [P,freqindex,timeindex]=one_half(data, wintype, winlen, step), 
%This function will calculate the 1 1/2 D Spectral surface. The 1 1/2 

%D surface characteristics are determined by the selection of window 

%type (wintype), window length (winlen) and the distance that the window 
%is moved through the data sequence (step). 

%The 1 % D surface (output matrix P) characteristics are determined by the selection 
Y%of window type (wintype), window length (winlen) and the distance that the 
Y%window is moved through the data sequence (step). 

%The program plots only the positive half of the spectral plane. The 
%outputs timeindex and freqindex are the appropriate plot indices for the 
%output time-frequency surface. 


% 

%The inputs are: 

Yodata: The input observations vector, for maximum effectiveness should 
% be of a length which is a power of 2, e.g. 64,128,512 
%wintype: 'O' Rectangular Window 

% 'l' Hamming Window 

%winlen: The desired width of the window, normally half of fhe input 
% length 

Ystep: Time step desired, can be ']' or a multiple of '2' 

% 

%The outputs are: 

%P The 1 % D time-frequency surface 


%timeindex The y axis index 
%freqindex The x axis index 

% ’ 
%See also ONESURF, ONELOFAR 


%Karen A. Hagerman 
%06 May 1992 


function [P,freqindex,timeindex}=one_half(data, wintype, winlen,step) 
[datarows,datacolumns]=size(data),; 
if datarows ~=1 
data=data.'; 
end 
siglen=length(data); 
if wintype==0 


win=ones(winlen-1,1); 
elseif wintype=1 
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win=hamming(winlen-1); 
end 


W=win(winlen/2:-1:1); 
x=[zeros(1,winlen) data zeros(1:winlen)].'; 
p=zeros(siglen/step, winlen); 


index=1; 

for n=winlen+1-step:siglent+winlen-step+1 
Xn=[abs(x(n))*2; abs(x(n))*2]; 
Xm=[conj(x(n:-1:n-(winlen/2-1))) x(n:n+(winlen/2-1))]; 
product=((Xm*Xn).*W).'; 
product=[product 0 conj(product(winlen/2:-1:2))]; 
p(index, :)=fftshift(abs(.5*ffRt(product))); 
index=index+1; 

end 

p=p(:,winlen/2+1 :winlen); 

[prow,pcolumn]=size(p); 


%Smoothing 

p_temp(1,:)=mean(p(1:3,:)); 

p_temp(2,:)=mean(p(1:4,:)); 

p_temp(prow-1,:)=mean(p(prow-3:prow,:)); 

p_temp(prow, :)}=mean(p(prow-2:prow,:)); 

for m=3:prow-2 
p_temp(m,:)}=mean(p(m-2:m+2,:)); 

end 


P=flipud(p_temp),; 


freqindex=[0:pcolumn-1]; 
timeindex=[1 :prow]; 
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ONESURF.M 


%function [Q,xindex, yindex]=onesurf{(sig, siglen, wintype, winlen, step); 
%This function will calculate an 1 1/2 D spectral surface made of 

%smaller 1 1/2 D surfaces. The smaller IPS surfaces are calculated for 
%data sequence pieces of length (siglen) of the total input data (sig). 
%Characteristics of the smaller surfaces are determined by the selection of 
%window type (wintype), window length (winlen) and the distance that the 
% window is moved through the data sequence (step). The surfaces are then 
%appended edge to edge to form the output Q surface. 


% 

%The inputs are: 

Ysig: The input observation sequence vector 

%sigien: The desired length of the discrete pieces of the input vector 
%wintype: 'O' Rectangular Window 

% ')' Hamming Window 

%winlen: The desired width of the window, normally half of the input 
% vector length 

%step: Time step desired, normally 'I' or a multiple of '2' 

% . 

%The outputs are: 

%P The ONESUFF time-frequency surface 


Yyindex The y axis index 
Yoxindex The x axis index 

% 

%See also ONE_HALF, ONELOFAR 


%Karen A. Hagerman 
%06 May 1992 


function [Q,xindex, yindex]=onesurf(data, siglen, wintype, winlen, step) 


if nargin==1 
siglen=input(Enter the desired length of the sequence '); 
wintype=input(ENTER "0" for RECT. WINDOW or "1" for HAMMING WINDOW::’); 
winlen=input((ENTER length of the window (must be an even number): '); 
step=input(Input desired step in time ‘); 
end 
[m,n]=size(data),; 
if m~=1 
data=data.'; 
end 
rows=floor(length(data)/siglen); 
data=[zeros(1,winlen) data zeros(1,winlen)]; 
len=length(data); 
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finish1=winlen/2-1; 
finish2=finish]+1; 


if wintype==0 
win=ones(winlen-1,1); 
elseif wintype=1 
win=hamming(winlen-1); 
end 
W=[win(winlen/2:-1:1)]'; 
Q=zeros(rows*(siglen/step),winlen/2), 
qindex=1:siglen/step:(rows*siglen/step); 
l=1; 
for k=1:siglen:len-siglen-2*winlen+1 
x=data(1,k:k+siglen+2*winlen-1); 
index=1; 
for n=winlen+1:step:siglen+winlen-step+ 1 
Xn=[abs(x(n))*2 abs(x(n))*2]; 
Xm=[conj(x(n:-1:n-(winlen/2-1))); x(n: oon 1): 
product=((Xn*Xm).*W); 
product=[product 0 conj(product(winlen/2:-1:2))]; 
P(index,:)=fftshift(real(.5*fit(product))); 
index=index+1; 
end 
Q(qindex(1):qindex(1)+(siglen/step)-1,:)=P(:,winlen/2+1-:winlen), 
l=I+1; 
end 
Q=flipud(abs(Q)); 
[qrow,qcolumn]=size(Q); 
xindex=[0:qcolumn- 1]; 
yindex=[1 :qrow]; 
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ONELOFAR.M 


%function [Q,freqindex, timeindex]=onelofar(sig, siglen, wintype, winlen, step), 
%This function will calculate a 'Lofar' display for a selected data sequence. 
%The total data sequence is first divided into equal length pieces of 
%length (siglen). An 1 1/2 D surface is then calculated for each piece. The 
%1 1/2 D surface characteristics are determined by the selection of window 
%type (wintype), window length (winlen) and the distance that the window 
%is moved through the data sequence pieces (step). The mean is then taken 
%of the 1 1/2 D surface and is placed sequentially in the Q matrix for display. 
%The Q matrix plots only the positive half of the spectral plane. The 
%outputs timeindex and freqindex can be used in plots to interpret the 
%results, 


% 
%The inputs are: 
%data: The input data string 


%siglen: The desired length of the discrete parts of hee sequence 
%wintype: 'O' Rectangular Window 


% '1!' Hamming Window 

%winlen: The desired width of the window, normally half of the siglen 
Ystep: Time step desired, normally '1' or a multiple of ‘2' 

%See also ONE_HALF, ONESURF 

% 

%The outputs are: 

%P The ONELOFAR time-frequency surface 


%timeindex The y axis index 
%freqindex The x axis index 
% 


%Karen A. Hagerman 
%06 May 1992 


function [Q, freqindex, timeindex]=onelofar(sig, siglen, wintype, winlen, step) 


if nargin==1 
siglen=input(Enter the desired length of the sequence '); 
wintype=input('ENTER "0" for RECT. WINDOW or "1" for HAMMING WINDOW:’); 
winlen=input(ENTER length of the window (must be an even number): '); 
step=input(‘Input desired step in time '); 

end 


{m,n]*size(sig); 

if m~—1 
sig=sig."; 

end 


59 








rows=floor(length(sig)/siglen); 
Q=zeros(rows, winlen/2); 
sig=[zeros(1,winlen) sig zeros(1,winlen)]; 
len=length(sig); 


if wintype==0 
win=ones(winlen-1,1); 
elseif wintype==1 
win=hamming(winlen-1); 

end 

W=win(winlen/2:-1:1).'; 


qindex=1; 
for k=1:siglen:len-siglen-2*winlen+1 
x=sig(1,k:k+siglen+(2*winlen)-1); 
index=1; 
for =winlen+1 :step:siglen+winlen-step+1 
Xn=[abs(x(1))*2 abs(x(1))*2]; 
Xm=[conj(x(1:-1:l-(winlen/2-1))); x: (winlen/2-1))]; 
product=((Xn* Xm).*W); 
product=[product 0 conj(product(winlen/2:-1:2))]; 
p(index, :)}=fitshift(abs(.5*fft(product))); 
index=index+1,; 
end 
if index ~=2 
Q(qindex, :)}=mean(p(:,winlen/2+1 :winlen)); 
gindex=qindex+1; 
else 
Q(qindex, :)=p(1,winlen/2+1 :winlen); 
qindex=qindex+1; 
end 
end 


Q=flipud(Q); 
[qrow,qcolumn]}=size(Q); 


freqindex=[0:qcolumn-1]; 
timeindex=[1:qrow]; 
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